#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include "calc_sobel.hpp"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"


__device__ inline void calcSobelOneInt(unsigned int * up, unsigned int * down, unsigned int *centre, float* dx_out, float * dy_out, float * dnorm_out)
{
	unsigned char left_char;
	unsigned char right_char;
	unsigned char up_char;
	unsigned char down_char;
//	unsigned char center_char;
	unsigned char left_up_char;
	unsigned char left_down_char;
	unsigned char right_up_char;
	unsigned char right_down_char;

	/*
	* for y orient:
	*  -1  -2  -1
	*   0   0   0
	*   1   2   1
	*/

	/*
	* for x orient:
	*  -1  0  1
	*  -2  0  2
	*  -1  0  1
	*/

	//1 calc the first byte
//	center_char = centre[1] & 0x00ff;
	left_char = (centre[0] & 0xff000000) >> 24;
	right_char = (centre[1] & 0x0000ff00) >> 8;
	up_char = up[1] & 0x00ff;
	down_char = down[1] & 0x00ff;

	left_up_char = (up[0] & 0xff000000) >> 24;
	left_down_char = (down[0] & 0xff000000) >> 24;
	right_up_char = (up[1] & 0x0000ff00) >> 8;
	right_down_char = (down[1] & 0x0000ff00) >> 8;

	dy_out[0] = float(left_down_char) + float((down_char << 1)) + float(right_down_char) - float(left_up_char) - float((up_char << 1)) - float(right_up_char);
	dx_out[0] = float(right_up_char) + float((right_char << 1)) + float(right_down_char) - float(left_up_char) - float((left_char << 1)) - float(left_down_char);
	dnorm_out[0] = dx_out[0] * dx_out[0] + dy_out[0] * dy_out[0];
	//2, calc the second byte

//	center_char = centre[1] & 0x0000ff00 >> 8;
	left_char = centre[1] & 0x00ff;
	right_char = (centre[1] & 0x00ff0000) >> 16;
	up_char = (up[1] & 0x0000ff00) >> 8;
	down_char = (down[1] & 0x0000ff00) >> 8;

	left_up_char = up[1] & 0x00ff;
	left_down_char = down[1] & 0x00ff;
	right_up_char = (up[1] & 0x00ff0000) >> 16;
	right_down_char = (down[1] & 0x00ff0000) >> 16;

	dy_out[1] = float(left_down_char) + float((down_char << 1)) + float(right_down_char) - float(left_up_char) - float((up_char << 1)) - float(right_up_char);
	//dx_out[1] = 58;
	dx_out[1] = float(right_up_char) + float((right_char << 1)) + float(right_down_char) - float(left_up_char) - float((left_char << 1)) - float(left_down_char);
	dnorm_out[1] = dx_out[1] * dx_out[1] + dy_out[1] * dy_out[1];

	//3, calc the third byte
//	center_char = centre[1] & 0x00ff0000 >> 16;
	left_char = (centre[1] >> 8) & 0x00ff;
	right_char = (centre[1] >> 24) & 0x00ff;
	up_char = (up[1] >> 16) & 0x00ff;
	down_char = (down[1] >> 16) & 0x00ff;

	left_up_char = (up[1] >> 8) & 0x00ff;
	left_down_char = (down[1] >> 8) & 0x00ff;
	right_up_char = (up[1] >> 24) & 0x00ff;
	right_down_char = (down[1] >> 24) & 0x00ff;
	dy_out[2] = float(left_down_char) + float((down_char << 1)) + float(right_down_char) - float(left_up_char) - float((up_char << 1)) - float(right_up_char);
	dx_out[2] = float(right_up_char) + float((right_char << 1)) + float(right_down_char) - float(left_up_char) - float((left_char << 1)) - float(left_down_char);
	dnorm_out[2] = dx_out[2] * dx_out[2] + dy_out[2] * dy_out[2];

	//4, calc the fourth byte
//	center_char = centre[1] & 0xff000000 >> 24;
	left_char = (centre[1] >> 16) & 0x00ff;
	right_char = centre[2] & 0x00ff;
	up_char = (up[1] >> 24) & 0x00ff;
	down_char = (down[1] >> 24) & 0x00ff;

	left_up_char = (up[1] >> 16) & 0x00ff;
	left_down_char = (down[1] >> 16) & 0x00ff;
	right_up_char = up[2] & 0x00ff;
	right_down_char = down[2] & 0x00ff;
	dy_out[3] = float(left_down_char) + float((down_char << 1)) + float(right_down_char) - float(left_up_char) - float((up_char << 1)) - float(right_up_char);
	dx_out[3] = float(right_up_char) + float((right_char << 1)) + float(right_down_char) - float(left_up_char) - float((left_char << 1)) - float(left_down_char);
	dnorm_out[3] = dx_out[3] * dx_out[3] + dy_out[3] * dy_out[3];

}

__global__ void sobel_64bit(unsigned int *d_in, float * dx_out,float * dy_out,float * dnorm_out, unsigned short width_64, unsigned short height, unsigned short with_8)
{
	//width_64: width of 8-bytes aligned;

	int x = blockDim.x*blockIdx.x + threadIdx.x;
	int y = blockDim.y*blockIdx.y + threadIdx.y;

	/*
	* in every thread, read 4 ints, that is 16 bytes.
	*
	* x is scaled between 0 -- width_64. increase x with 1, equalize index for d_in increase with 2.
	*/

	if ((x >= width_64) || (y >= height))
	{
		return;
	}
	//d_in is type of int, so x should multiply with 2.
	int index = y*(width_64<<1) + (x<<1);

	if (x < (width_64 - 2))
	{
		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height-1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height-1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height-1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height-1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		if (x == 0)
		{
			tempUp[0] = (tempUpInt0 & 0x00ff00) << 16;
			tempUp[1] = tempUpInt0;
			tempUp[2] = tempUpInt1;

			tempCentre[0] = (tempInt0 & 0x00ff00) << 16;
			tempCentre[1] = tempInt0;
			tempCentre[2] = tempInt1;

			tempDown[0] = (tempDownInt0 & 0x00ff00) << 16;
			tempDown[1] = tempDownInt0;
			tempDown[2] = tempDownInt1;

			calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2), dy_out + (index << 2), dnorm_out + (index << 2));
		}

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		tempUp[0] = tempUpInt1;
		tempUp[1] = tempUpInt2;
		tempUp[2] = tempUpInt3;

		tempCentre[0] = tempInt1;
		tempCentre[1] = tempInt2;
		tempCentre[2] = tempInt3;

		tempDown[0] = tempDownInt1;
		tempDown[1] = tempDownInt2;
		tempDown[2] = tempDownInt3;

		calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);
	}
	else if (x == (width_64 - 2))
	{

		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		unsigned short remain = with_8 - (x+1) * 8;

		if(remain==0)
		{
             return;
		}

		if(remain > 4)
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = tempUpInt3;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = tempInt3;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = tempDownInt3;
		}
		else if(remain == 4)
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = (tempUpInt2 >> 16)&0x00ff;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = (tempInt2 >> 16)&0x00ff;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = (tempDownInt2 >> 16)&0x00ff;
		}
		else if(remain == 1)
		{

			tempUp[0] = tempUpInt1;
			tempUp[1] = (tempUpInt2&0xff)|((tempUpInt1>>24)<<8);
			tempUp[2] = 0;

			tempCentre[0] = tempInt1;
			tempCentre[1] = (tempInt2&0xff)|((tempInt1>>24)<<8);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt1;
			tempDown[1] = (tempDownInt2&0xff)|((tempDownInt1>>24)<<8);
			tempDown[2] = 0;
		}
		else
		{
			unsigned int mask = 0xffffffff>>((4-remain)*8);
			unsigned int mask1 = 0x00ff <<((remain-2)*8);
			tempUp[0] = tempUpInt1;
			tempUp[1] = (tempUpInt2&mask)|((tempUpInt2&mask1)<<16);
			tempUp[2] = 0;

			tempCentre[0] = tempInt1;
			tempCentre[1] = (tempInt2&mask)|((tempInt2&mask1)<<16);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt1;
			tempDown[1] = (tempDownInt2&mask)|((tempDownInt2&mask1)<<16);
			tempDown[2] = 0;
		}


		calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);

	}
	else if (x == (width_64 - 1))
	{
		unsigned short remain = with_8 - x * 8;
		/*if remain < 4, do nothing*/
		if (remain > 4)
		{
			unsigned int tempInt0 = d_in[index];
			unsigned int tempInt1 = d_in[index + 1];

			unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
			unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];

			unsigned int tempDownInt0 = (y == height-1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
			unsigned int tempDownInt1 = (y == height-1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];

			unsigned int tempUp[3];
			unsigned int tempCentre[3];
			unsigned int tempDown[3];

			if(remain == 5)
			{
				tempUp[0] = tempUpInt0;
				tempUp[1] = (tempUpInt1&0x00ff)|((tempUpInt0>>24)<<8);
				tempUp[2] = 0;

				tempCentre[0] = tempInt0;
				tempCentre[1] = (tempInt1&0x00ff)|((tempInt0>>24)<<8);
				tempCentre[2] = 0;

				tempDown[0] = tempDownInt0;
				tempDown[1] = (tempDownInt1&0x00ff)|((tempDownInt0>>24)<<8);
				tempDown[2] = 0;
			}
			else if(remain == 8)
			{
				tempUp[0] = tempUpInt0;
				tempUp[1] = tempUpInt1;
				tempUp[2] = (tempUpInt1>>16)&0x00ff;

				tempCentre[0] = tempInt0;
				tempCentre[1] = tempInt1;
				tempCentre[2] = (tempInt1>>16)&0x00ff;

				tempDown[0] = tempDownInt0;
				tempDown[1] = tempDownInt1;
				tempDown[2] = (tempDownInt1>>16)&0x00ff;
			}
			else
			{
				unsigned int mask = 0xffffffff>>((8-remain)*8);
				unsigned int mask1 = 0x00ff <<((remain-6)*8);

	            tempUp[0] = tempUpInt0;
				tempUp[1] = (tempUpInt1&mask)|((tempUpInt1&mask1)<<16);
				tempUp[2] = 0;

				tempCentre[0] = tempInt0;
				tempCentre[1] = (tempInt1&mask)|((tempInt1&mask1)<<16);
				tempCentre[2] = 0;

				tempDown[0] = tempDownInt0;
				tempDown[1] = (tempDownInt1&mask)|((tempDownInt1&mask1)<<16);
				tempDown[2] = 0;
			}

			calcSobelOneInt(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);
		}
	}
}


int calc_sobel_cuda(uchar * dev_imageBuffIn, unsigned int width, unsigned int height, uchar * devBuff, int pitch, float * host_dx, float * host_dy, float *host_norm)
{
	if ((dev_imageBuffIn == NULL) || (host_dx == NULL) || (host_dy == NULL) || (host_norm == NULL))
	{
		return -1;
	}

	dim3 blockSize2(TX, TY);
	int blockNumY = (height + TY - 1) / TY;
	int blockNumX2 = ((width + 7) / 8 + TX - 1) / TX;
	dim3 gridSize2(blockNumX2, blockNumY);

	float * dxOut = (float*)devBuff;

	float * dyOut = (float*)(devBuff + sizeof(float)*pitch*height);

	float * dnormOut = (float*)(devBuff + ((sizeof(float)*pitch*height) << 1));

    sobel_64bit << < gridSize2, blockSize2 >> > ((unsigned int*)dev_imageBuffIn, dxOut, dyOut, dnormOut, (width + 7) >> 3, height, width);

	int sizeFloat = sizeof(float);

	if (width == pitch)
	{
		cudaMemcpy(host_dx, dxOut, height*width*sizeFloat, cudaMemcpyDeviceToHost);
		cudaMemcpy(host_dy, dyOut, height*width*sizeFloat, cudaMemcpyDeviceToHost);
		cudaMemcpy(host_norm, dnormOut, height*width*sizeFloat, cudaMemcpyDeviceToHost);
	}
	else
	{
		cudaMemcpy2D(host_dx, width*sizeFloat, dxOut, pitch*sizeFloat, width*sizeFloat, height, cudaMemcpyDeviceToHost);
		cudaMemcpy2D(host_dy, width*sizeFloat, dyOut, pitch*sizeFloat, width*sizeFloat, height, cudaMemcpyDeviceToHost);
		cudaMemcpy2D(host_norm, width*sizeFloat, dnormOut, pitch*sizeFloat, width*sizeFloat, height, cudaMemcpyDeviceToHost);
	}
	return 0;
}

__device__ inline void calcSobelOneInt_short(unsigned int * up, unsigned int * down, unsigned int *centre, short * dx_out, short * dy_out, int * dnorm_out)
{
	unsigned char left_char;
	unsigned char right_char;
	unsigned char up_char;
	unsigned char down_char;
	//	unsigned char center_char;
	unsigned char left_up_char;
	unsigned char left_down_char;
	unsigned char right_up_char;
	unsigned char right_down_char;

	/*
	* for y orient:
	*  -1  -2  -1
	*   0   0   0
	*   1   2   1
	*/

	/*
	* for x orient:
	*  -1  0  1
	*  -2  0  2
	*  -1  0  1
	*/

	//1 calc the first byte
	//	center_char = centre[1] & 0x00ff;
	left_char = (centre[0] & 0xff000000) >> 24;
	right_char = (centre[1] & 0x0000ff00) >> 8;
	up_char = up[1] & 0x00ff;
	down_char = down[1] & 0x00ff;

	left_up_char = (up[0] & 0xff000000) >> 24;
	left_down_char = (down[0] & 0xff000000) >> 24;
	right_up_char = (up[1] & 0x0000ff00) >> 8;
	right_down_char = (down[1] & 0x0000ff00) >> 8;

	dy_out[0] = short(left_down_char) + short((down_char << 1)) + short(right_down_char) - short(left_up_char) - short((up_char << 1)) - short(right_up_char);
	dx_out[0] = short(right_up_char) + short((right_char << 1)) + short(right_down_char) - short(left_up_char) - short((left_char << 1)) - short(left_down_char);
	dnorm_out[0] = int(dx_out[0] * dx_out[0]) + int(dy_out[0] * dy_out[0]);
	//2, calc the second byte

	//	center_char = centre[1] & 0x0000ff00 >> 8;
	left_char = centre[1] & 0x00ff;
	right_char = (centre[1] & 0x00ff0000) >> 16;
	up_char = (up[1] & 0x0000ff00) >> 8;
	down_char = (down[1] & 0x0000ff00) >> 8;

	left_up_char = up[1] & 0x00ff;
	left_down_char = down[1] & 0x00ff;
	right_up_char = (up[1] & 0x00ff0000) >> 16;
	right_down_char = (down[1] & 0x00ff0000) >> 16;

	dy_out[1] = short(left_down_char) + short((down_char << 1)) + short(right_down_char) - short(left_up_char) - short((up_char << 1)) - short(right_up_char);
	//dx_out[1] = 58;
	dx_out[1] = short(right_up_char) + short((right_char << 1)) + short(right_down_char) - short(left_up_char) - short((left_char << 1)) - short(left_down_char);
	dnorm_out[1] = int(dx_out[1] * dx_out[1]) + int(dy_out[1] * dy_out[1]);

	//3, calc the third byte
	//	center_char = centre[1] & 0x00ff0000 >> 16;
	left_char = (centre[1] >> 8) & 0x00ff;
	right_char = (centre[1] >> 24) & 0x00ff;
	up_char = (up[1] >> 16) & 0x00ff;
	down_char = (down[1] >> 16) & 0x00ff;

	left_up_char = (up[1] >> 8) & 0x00ff;
	left_down_char = (down[1] >> 8) & 0x00ff;
	right_up_char = (up[1] >> 24) & 0x00ff;
	right_down_char = (down[1] >> 24) & 0x00ff;
	dy_out[2] = short(left_down_char) + short((down_char << 1)) + short(right_down_char) - short(left_up_char) - short((up_char << 1)) - short(right_up_char);
	dx_out[2] = short(right_up_char) + short((right_char << 1)) + short(right_down_char) - short(left_up_char) - short((left_char << 1)) - short(left_down_char);
	dnorm_out[2] = int(dx_out[2] * dx_out[2]) + int(dy_out[2] * dy_out[2]);

	//4, calc the fourth byte
	//	center_char = centre[1] & 0xff000000 >> 24;
	left_char = (centre[1] >> 16) & 0x00ff;
	right_char = centre[2] & 0x00ff;
	up_char = (up[1] >> 24) & 0x00ff;
	down_char = (down[1] >> 24) & 0x00ff;

	left_up_char = (up[1] >> 16) & 0x00ff;
	left_down_char = (down[1] >> 16) & 0x00ff;
	right_up_char = up[2] & 0x00ff;
	right_down_char = down[2] & 0x00ff;
	dy_out[3] = short(left_down_char) + short((down_char << 1)) + short(right_down_char) - short(left_up_char) - short((up_char << 1)) - short(right_up_char);
	dx_out[3] = short(right_up_char) + short((right_char << 1)) + short(right_down_char) - short(left_up_char) - short((left_char << 1)) - short(left_down_char);
	dnorm_out[3] = int(dx_out[3] * dx_out[3]) + int(dy_out[3] * dy_out[3]);

}

__global__ void sobel_64bit_short(unsigned int *d_in, short * dx_out, short * dy_out, int * dnorm_out, unsigned short width_64, unsigned short height, unsigned short with_8)
{
	//width_64: width of 8-bytes aligned;

	int x = blockDim.x*blockIdx.x + threadIdx.x;
	int y = blockDim.y*blockIdx.y + threadIdx.y;

	/*
	* in every thread, read 4 ints, that is 16 bytes.
	*
	* x is scaled between 0 -- width_64. increase x with 1, equalize index for d_in increase with 2.
	*/

	if ((x >= width_64) || (y >= height))
	{
		return;
	}
	//d_in is type of int, so x should multiply with 2.
	int index = y*(width_64 << 1) + (x << 1);

	if (x < (width_64 - 2))
	{
		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		if (x == 0)
		{
			tempUp[0] = (tempUpInt0 & 0x00ff00) << 16;
			tempUp[1] = tempUpInt0;
			tempUp[2] = tempUpInt1;

			tempCentre[0] = (tempInt0 & 0x00ff00) << 16;
			tempCentre[1] = tempInt0;
			tempCentre[2] = tempInt1;

			tempDown[0] = (tempDownInt0 & 0x00ff00) << 16;
			tempDown[1] = tempDownInt0;
			tempDown[2] = tempDownInt1;

			calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2), dy_out + (index << 2), dnorm_out + (index << 2));
		}

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		tempUp[0] = tempUpInt1;
		tempUp[1] = tempUpInt2;
		tempUp[2] = tempUpInt3;

		tempCentre[0] = tempInt1;
		tempCentre[1] = tempInt2;
		tempCentre[2] = tempInt3;

		tempDown[0] = tempDownInt1;
		tempDown[1] = tempDownInt2;
		tempDown[2] = tempDownInt3;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);
	}
	else if (x == (width_64 - 2))
	{

		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		unsigned short remain = with_8 - (x + 1) * 8;

		if (remain == 0)
		{
			return;
		}

		if (remain > 4)
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = tempUpInt3;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = tempInt3;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = tempDownInt3;
		}
		else if (remain == 4)
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = (tempUpInt2 >> 16) & 0x00ff;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = (tempInt2 >> 16) & 0x00ff;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = (tempDownInt2 >> 16) & 0x00ff;
		}
		else if (remain == 1)
		{

			tempUp[0] = tempUpInt1;
			tempUp[1] = (tempUpInt2 & 0xff) | ((tempUpInt1 >> 24) << 8);
			tempUp[2] = 0;

			tempCentre[0] = tempInt1;
			tempCentre[1] = (tempInt2 & 0xff) | ((tempInt1 >> 24) << 8);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt1;
			tempDown[1] = (tempDownInt2 & 0xff) | ((tempDownInt1 >> 24) << 8);
			tempDown[2] = 0;
		}
		else
		{
			unsigned int mask = 0xffffffff >> ((4 - remain) * 8);
			unsigned int mask1 = 0x00ff << ((remain - 2) * 8);
			tempUp[0] = tempUpInt1;
			tempUp[1] = (tempUpInt2&mask) | ((tempUpInt2&mask1) << 16);
			tempUp[2] = 0;

			tempCentre[0] = tempInt1;
			tempCentre[1] = (tempInt2&mask) | ((tempInt2&mask1) << 16);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt1;
			tempDown[1] = (tempDownInt2&mask) | ((tempDownInt2&mask1) << 16);
			tempDown[2] = 0;
		}


		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);

	}
	else if (x == (width_64 - 1))
	{
		unsigned short remain = with_8 - x * 8;
		/*if remain < 4, do nothing*/
		if (remain > 4)
		{
			unsigned int tempInt0 = d_in[index];
			unsigned int tempInt1 = d_in[index + 1];

			unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
			unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];

			unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
			unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];

			unsigned int tempUp[3];
			unsigned int tempCentre[3];
			unsigned int tempDown[3];

			if (remain == 5)
			{
				tempUp[0] = tempUpInt0;
				tempUp[1] = (tempUpInt1 & 0x00ff) | ((tempUpInt0 >> 24) << 8);
				tempUp[2] = 0;

				tempCentre[0] = tempInt0;
				tempCentre[1] = (tempInt1 & 0x00ff) | ((tempInt0 >> 24) << 8);
				tempCentre[2] = 0;

				tempDown[0] = tempDownInt0;
				tempDown[1] = (tempDownInt1 & 0x00ff) | ((tempDownInt0 >> 24) << 8);
				tempDown[2] = 0;
			}
			else if (remain == 8)
			{
				tempUp[0] = tempUpInt0;
				tempUp[1] = tempUpInt1;
				tempUp[2] = (tempUpInt1 >> 16) & 0x00ff;

				tempCentre[0] = tempInt0;
				tempCentre[1] = tempInt1;
				tempCentre[2] = (tempInt1 >> 16) & 0x00ff;

				tempDown[0] = tempDownInt0;
				tempDown[1] = tempDownInt1;
				tempDown[2] = (tempDownInt1 >> 16) & 0x00ff;
			}
			else
			{
				unsigned int mask = 0xffffffff >> ((8 - remain) * 8);
				unsigned int mask1 = 0x00ff << ((remain - 6) * 8);

				tempUp[0] = tempUpInt0;
				tempUp[1] = (tempUpInt1&mask) | ((tempUpInt1&mask1) << 16);
				tempUp[2] = 0;

				tempCentre[0] = tempInt0;
				tempCentre[1] = (tempInt1&mask) | ((tempInt1&mask1) << 16);
				tempCentre[2] = 0;

				tempDown[0] = tempDownInt0;
				tempDown[1] = (tempDownInt1&mask) | ((tempDownInt1&mask1) << 16);
				tempDown[2] = 0;
			}

			calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);
		}
	}
}


__global__ void sobel_64bit_short_v1(unsigned int *d_in, short * dx_out, short * dy_out, int * dnorm_out, unsigned short width_64, unsigned short height, unsigned short with_8)
{
	//width_64: width of 8-bytes aligned;

	int x = blockDim.x*blockIdx.x + threadIdx.x;
	int y = blockDim.y*blockIdx.y + threadIdx.y;

	/*
	* in every thread, read 4 ints, that is 16 bytes.
	*
	* x is scaled between 0 -- width_64. increase x with 1, equalize index for d_in increase with 2.
	*/

	if ((x >= width_64) || (y >= height))
	{
		return;
	}
	//d_in is type of int, so x should multiply with 2.
	int index = y*(width_64 << 1) + (x << 1);

	if (x < (width_64 - 2))
	{
		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		if (x == 0)
		{
			tempUp[0] = (tempUpInt0 & 0x00ff00) << 16;
			tempUp[1] = tempUpInt0;
			tempUp[2] = tempUpInt1;

			tempCentre[0] = (tempInt0 & 0x00ff00) << 16;
			tempCentre[1] = tempInt0;
			tempCentre[2] = tempInt1;

			tempDown[0] = (tempDownInt0 & 0x00ff00) << 16;
			tempDown[1] = tempDownInt0;
			tempDown[2] = tempDownInt1;

			calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2), dy_out + (index << 2), dnorm_out + (index << 2));
		}

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		tempUp[0] = tempUpInt1;
		tempUp[1] = tempUpInt2;
		tempUp[2] = tempUpInt3;

		tempCentre[0] = tempInt1;
		tempCentre[1] = tempInt2;
		tempCentre[2] = tempInt3;

		tempDown[0] = tempDownInt1;
		tempDown[1] = tempDownInt2;
		tempDown[2] = tempDownInt3;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);
	}
	else if (x == (width_64 - 2))
	{

		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		unsigned short remain = with_8 - (x + 1) * 8;

		if (remain == 0)
		{
			return;
		}

		if (remain > 4)
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = tempUpInt3;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = tempInt3;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = tempDownInt3;
		}
		else
		{
			tempUp[0] = tempUpInt1;
			tempUp[1] = tempUpInt2;
			tempUp[2] = 0;

			tempCentre[0] = tempInt1;
			tempCentre[1] = tempInt2;
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt1;
			tempDown[1] = tempDownInt2;
			tempDown[2] = 0;
		}

		/*if (remain == 4)
		{
		tempUp[0] = tempUpInt1;
		tempUp[1] = tempUpInt2;
		tempUp[2] = (tempUpInt2 >> 16) & 0x00ff;

		tempCentre[0] = tempInt1;
		tempCentre[1] = tempInt2;
		tempCentre[2] = (tempInt2 >> 16) & 0x00ff;

		tempDown[0] = tempDownInt1;
		tempDown[1] = tempDownInt2;
		tempDown[2] = (tempDownInt2 >> 16) & 0x00ff;
		}
		else if (remain == 1)
		{

		tempUp[0] = tempUpInt1;
		tempUp[1] = (tempUpInt2 & 0xff) | ((tempUpInt1 >> 24) << 8);
		tempUp[2] = 0;

		tempCentre[0] = tempInt1;
		tempCentre[1] = (tempInt2 & 0xff) | ((tempInt1 >> 24) << 8);
		tempCentre[2] = 0;

		tempDown[0] = tempDownInt1;
		tempDown[1] = (tempDownInt2 & 0xff) | ((tempDownInt1 >> 24) << 8);
		tempDown[2] = 0;
		}
		else
		{
		unsigned int mask = 0xffffffff >> ((4 - remain) * 8);
		unsigned int mask1 = 0x00ff << ((remain - 2) * 8);
		tempUp[0] = tempUpInt1;
		tempUp[1] = (tempUpInt2&mask) | ((tempUpInt2&mask1) << 16);
		tempUp[2] = 0;

		tempCentre[0] = tempInt1;
		tempCentre[1] = (tempInt2&mask) | ((tempInt2&mask1) << 16);
		tempCentre[2] = 0;

		tempDown[0] = tempDownInt1;
		tempDown[1] = (tempDownInt2&mask) | ((tempDownInt2&mask1) << 16);
		tempDown[2] = 0;
		}*/


		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);

	}
	else if (x == (width_64 - 1))
	{
		unsigned short remain = with_8 - x * 8;
		/*if remain < 4, do nothing*/
		if (remain > 4)
		{
			unsigned int tempInt0 = d_in[index];
			unsigned int tempInt1 = d_in[index + 1];

			unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
			unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];

			unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
			unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];

			unsigned int tempUp[3];
			unsigned int tempCentre[3];
			unsigned int tempDown[3];


			tempUp[0] = tempUpInt0;
			tempUp[1] = tempUpInt1;
			tempUp[2] = 0;

			tempCentre[0] = tempInt0;
			tempCentre[1] = tempInt1;
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt0;
			tempDown[1] = tempDownInt1;
			tempDown[2] = 0;

			/*if (remain == 5)
			{
			tempUp[0] = tempUpInt0;
			tempUp[1] = (tempUpInt1 & 0x00ff) | ((tempUpInt0 >> 24) << 8);
			tempUp[2] = 0;

			tempCentre[0] = tempInt0;
			tempCentre[1] = (tempInt1 & 0x00ff) | ((tempInt0 >> 24) << 8);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt0;
			tempDown[1] = (tempDownInt1 & 0x00ff) | ((tempDownInt0 >> 24) << 8);
			tempDown[2] = 0;
			}
			else if (remain == 8)
			{
			tempUp[0] = tempUpInt0;
			tempUp[1] = tempUpInt1;
			tempUp[2] = (tempUpInt1 >> 16) & 0x00ff;

			tempCentre[0] = tempInt0;
			tempCentre[1] = tempInt1;
			tempCentre[2] = (tempInt1 >> 16) & 0x00ff;

			tempDown[0] = tempDownInt0;
			tempDown[1] = tempDownInt1;
			tempDown[2] = (tempDownInt1 >> 16) & 0x00ff;
			}
			else
			{
			unsigned int mask = 0xffffffff >> ((8 - remain) * 8);
			unsigned int mask1 = 0x00ff << ((remain - 6) * 8);

			tempUp[0] = tempUpInt0;
			tempUp[1] = (tempUpInt1&mask) | ((tempUpInt1&mask1) << 16);
			tempUp[2] = 0;

			tempCentre[0] = tempInt0;
			tempCentre[1] = (tempInt1&mask) | ((tempInt1&mask1) << 16);
			tempCentre[2] = 0;

			tempDown[0] = tempDownInt0;
			tempDown[1] = (tempDownInt1&mask) | ((tempDownInt1&mask1) << 16);
			tempDown[2] = 0;
			}*/

			calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);
		}
	}
}


int calc_sobel_cuda(uchar * dev_imageBuffIn, unsigned int width, unsigned int height, uchar * devBuff, int pitch, short * host_dx, short * host_dy, int *host_norm)
{
	//uchar * d_in;

	if ((dev_imageBuffIn == NULL) || (host_dx == NULL) || (host_dy == NULL) || (host_norm == NULL))
	{
		return -1;
	}

	//d_in = devBuff;

	dim3 blockSize2(TX, TY);
	int blockNumY = (height + TY - 1) / TY;
	int blockNumX2 = ((width + 7) / 8 + TX - 1) / TX;
	dim3 gridSize2(blockNumX2, blockNumY);

	short * dxOut = (short*)devBuff;

	short * dyOut = (short*)(dxOut + pitch*height);

	int * dnormOut = (int*)(dyOut + pitch*height);

	//printf("gridSize2.x=%d,gridSize2.y=%d\n", gridSize2.x, gridSize2.y);

	sobel_64bit_short << < gridSize2, blockSize2 >> > ((unsigned int*)dev_imageBuffIn, dxOut, dyOut, dnormOut, (width + 7) >> 3, height, width);

	//	int sizeFloat = sizeof(float);
	int sizeInt = sizeof(int);
	int sizeShort = sizeof(short);

	if (width == pitch)
	{
		cudaMemcpy(host_dx, dxOut, height*width*sizeShort, cudaMemcpyDeviceToHost);
		cudaMemcpy(host_dy, dyOut, height*width*sizeShort, cudaMemcpyDeviceToHost);
		cudaMemcpy(host_norm, dnormOut, height*width*sizeInt, cudaMemcpyDeviceToHost);
	}
	else
	{
		cudaMemcpy2D(host_dx, width*sizeShort, dxOut, pitch*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);
		cudaMemcpy2D(host_dy, width*sizeShort, dyOut, pitch*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);
		cudaMemcpy2D(host_norm, width*sizeInt, dnormOut, pitch*sizeInt, width*sizeInt, height, cudaMemcpyDeviceToHost);
	}
	return 0;
}


__global__ void sobel_64bit_short_v2(unsigned int *d_in, short * dx_out, short * dy_out, int * dnorm_out, unsigned short width_64, unsigned short height, unsigned short with_8)
{
	//width_64: width of 8-bytes aligned;

	int x = blockDim.x*blockIdx.x + threadIdx.x;
	int y = blockDim.y*blockIdx.y + threadIdx.y;

	/*
	* in every thread, read 4 ints, that is 16 bytes.
	*
	* x is scaled between 0 -- width_64. increase x with 1, equalize index for d_in increase with 2.
	*/

	if ((x >= width_64) || (y >= height))
	{
		return;
	}
	//d_in is type of int, so x should multiply with 2.
	int index = y*(width_64 << 1) + (x << 1);

	if (x < (width_64 - 1))
	{
		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = d_in[index + 2];
		unsigned int tempInt3 = d_in[index + 3];

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = (y == 0) ? d_in[index + (width_64 << 1) + 2] : d_in[index - (width_64 << 1) + 2];
		unsigned int tempUpInt3 = (y == 0) ? d_in[index + (width_64 << 1) + 3] : d_in[index - (width_64 << 1) + 3];

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = (y == height - 1) ? tempUpInt2 : d_in[index + (width_64 << 1) + 2];
		unsigned int tempDownInt3 = (y == height - 1) ? tempUpInt3 : d_in[index + (width_64 << 1) + 3];

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		if (x == 0)
		{
			tempUp[0] = 0;
			tempUp[1] = tempUpInt0;
			tempUp[2] = tempUpInt1;

			tempCentre[0] = 0;
			tempCentre[1] = tempInt0;
			tempCentre[2] = tempInt1;

			tempDown[0] = 0;
			tempDown[1] = tempDownInt0;
			tempDown[2] = tempDownInt1;

			calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2), dy_out + (index << 2), dnorm_out + (index << 2));
		}

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

		tempUp[0] = tempUpInt1;
		tempUp[1] = tempUpInt2;
		tempUp[2] = tempUpInt3;

		tempCentre[0] = tempInt1;
		tempCentre[1] = tempInt2;
		tempCentre[2] = tempInt3;

		tempDown[0] = tempDownInt1;
		tempDown[1] = tempDownInt2;
		tempDown[2] = tempDownInt3;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 8, dy_out + (index << 2) + 8, dnorm_out + (index << 2) + 8);
	}
	else if (x == (width_64 - 1))
	{

		unsigned int tempInt0 = d_in[index];
		unsigned int tempInt1 = d_in[index + 1];
		unsigned int tempInt2 = 0;

		unsigned int tempUpInt0 = (y == 0) ? d_in[index + (width_64 << 1)] : d_in[index - (width_64 << 1)];
		unsigned int tempUpInt1 = (y == 0) ? d_in[index + (width_64 << 1) + 1] : d_in[index - (width_64 << 1) + 1];
		unsigned int tempUpInt2 = 0;

		unsigned int tempDownInt0 = (y == height - 1) ? tempUpInt0 : d_in[index + (width_64 << 1)];
		unsigned int tempDownInt1 = (y == height - 1) ? tempUpInt1 : d_in[index + (width_64 << 1) + 1];
		unsigned int tempDownInt2 = 0;

		unsigned int tempUp[3];
		unsigned int tempCentre[3];
		unsigned int tempDown[3];

		tempUp[0] = tempUpInt0;
		tempUp[1] = tempUpInt1;
		tempUp[2] = tempUpInt2;

		tempCentre[0] = tempInt0;
		tempCentre[1] = tempInt1;
		tempCentre[2] = tempInt2;

		tempDown[0] = tempDownInt0;
		tempDown[1] = tempDownInt1;
		tempDown[2] = tempDownInt2;

		calcSobelOneInt_short(tempUp, tempDown, tempCentre, dx_out + (index << 2) + 4, dy_out + (index << 2) + 4, dnorm_out + (index << 2) + 4);

	}

}


int calc_sobel_cuda_v1(uchar * dev_imageBuffIn, unsigned int width, unsigned int height, uchar * devBuff, short * host_dx, short * host_dy, int *host_norm)
{
	//uchar * d_in;

	if ((dev_imageBuffIn == NULL) || (host_dx == NULL) || (host_dy == NULL) || (host_norm == NULL))
	{
		return -1;
	}
	uchar * devImgExpanded = NULL;

	expand_border_dev(dev_imageBuffIn, width, height, devBuff, &devImgExpanded);

	//d_in = devBuff;

	int widthEx = width + 2;
	int heightEx = height + 2;

	int pitch1 = ((widthEx + 7) >> 3) << 3;

	dim3 blockSize2(TX, TY);
	int blockNumY = (heightEx + TY - 1) / TY;
	int blockNumX2 = ((widthEx + 7) / 8 + TX - 1) / TX;
	dim3 gridSize2(blockNumX2, blockNumY);

	uchar * devImgExpandedPitched = devImgExpanded + heightEx * pitch1;

	cudaMemcpy2D(devImgExpandedPitched, pitch1, devImgExpanded, widthEx, widthEx, heightEx,cudaMemcpyDeviceToDevice);

	short * dxOut = (short*)(devImgExpandedPitched + pitch1*heightEx);

	short * dyOut = (short*)(dxOut + pitch1*heightEx);

	int * dnormOut = (int*)(dyOut + pitch1*heightEx);

	//printf("gridSize2.x=%d,gridSize2.y=%d\n", gridSize2.x, gridSize2.y);

	sobel_64bit_short_v2 << < gridSize2, blockSize2 >> > ((unsigned int*)devImgExpandedPitched, dxOut, dyOut, dnormOut, (widthEx + 7) >> 3, heightEx, widthEx);

	int sizeFloat = sizeof(float);
	int sizeInt = sizeof(int);
	int sizeShort = sizeof(short);

	/*if (width == pitch)
	{
	cudaMemcpy(host_dx, dxOut, height*width*sizeShort, cudaMemcpyDeviceToHost);
	cudaMemcpy(host_dy, dyOut, height*width*sizeShort, cudaMemcpyDeviceToHost);
	cudaMemcpy(host_norm, dnormOut, height*width*sizeInt, cudaMemcpyDeviceToHost);
	}
	else
	{
	cudaMemcpy2D(host_dx, width*sizeShort, dxOut, pitch*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);
	cudaMemcpy2D(host_dy, width*sizeShort, dyOut, pitch*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);
	cudaMemcpy2D(host_norm, width*sizeInt, dnormOut, pitch*sizeInt, width*sizeInt, height, cudaMemcpyDeviceToHost);
	}*/

	cudaMemcpy2D(host_dx, width*sizeShort, dxOut + pitch1 + 1, pitch1*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);

	//cudaMemcpy(host_dx, dxOut, pitch1*sizeShort*4, cudaMemcpyDeviceToHost);
	
	//cudaMemcpy2D(host_dy, widthEx, devImgExpandedPitched, pitch1, widthEx, 4, cudaMemcpyDeviceToHost);
	cudaMemcpy2D(host_dy, width*sizeShort, dyOut + pitch1 + 1, pitch1*sizeShort, width*sizeShort, height, cudaMemcpyDeviceToHost);
	cudaMemcpy2D(host_norm, width*sizeInt, dnormOut + pitch1 + 1, pitch1*sizeInt, width*sizeInt, height, cudaMemcpyDeviceToHost);
	return 0;
}



__global__ void expand_border_kernel(uchar* devInOut, unsigned int width, unsigned int height)
{
	if ((blockIdx.x >= height) || (threadIdx.x >= 2))
	{
		return;
	}

	int index = blockIdx.x * width + threadIdx.x * (width - 1);

	if (threadIdx.x == 0)
	{
		devInOut[index] = devInOut[index + 2];
	}
	else
	{
		devInOut[index] = devInOut[index - 2];
	}

}


__global__ void expand_border_kernel_v1(uchar* devInOut, unsigned int width, unsigned int height)
{
	int indexLine = (blockIdx.x << 5) + threadIdx.x;

	if ((indexLine >= height) || (threadIdx.y >= 2))
	{
		return;
	}

	int indexLineHead = indexLine * width;
	int indexDst = (threadIdx.y == 0) ? indexLineHead : indexLineHead + width - 1;
	int indexSrc = (threadIdx.y == 0) ? indexDst + 2 : indexDst - 2;

	devInOut[indexDst] = devInOut[indexSrc];
}


int expand_border_dev(uchar * dev_imgIn, unsigned int widthSrc, unsigned int heightSrc, uchar * devBuff, uchar** hostImgOut)
{
	dim3 blockSize(32, 32);
	dim3 gridSize((widthSrc + 31 / 32), (heightSrc + 31) / 32);

	unsigned int widthDst = widthSrc + 2;
	unsigned int heightDst = heightSrc + 2;

	cudaMemcpy(devBuff + 1, dev_imgIn + widthSrc, widthSrc, cudaMemcpyDeviceToDevice);
	cudaMemcpy(devBuff + 1 + widthDst * (heightSrc + 1), dev_imgIn + widthSrc*(heightSrc - 2), widthSrc, cudaMemcpyDeviceToDevice);

	cudaMemcpy2D(devBuff + 1 + widthDst, widthDst, dev_imgIn, widthSrc, widthSrc, heightSrc, cudaMemcpyDeviceToDevice);
	//expand_border_kernel << <heightDst, 2 >> >(devBuff, widthDst, heightDst);
	dim3 blockSize1(32, 2);
	int gridX = (heightDst + 31) / 32;
	expand_border_kernel_v1 << <gridX, blockSize1 >> >(devBuff, widthDst, heightDst);
	
	*hostImgOut = devBuff;

	return 0;
}

